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ABSTRACT 

We present a statistical study of the effects induced by substructures on the deflection 
potential of dark matter halos in the strong lensing regime. This investigation is based 
on the pertubative solution around the Einstein radius (Alard 2007) in which all the 
information on the deflection potential is specified by only a pair of one-dimensional 
functions on this ring. 

Using direct comparison with ray-tracing solutions, we found that the iso-contours 
of lensed images predicted by the pertubative solution is reproduced with a mean 
error on their radial extension of less than 1% — in units of the Einstein radius, for 
reasonable substructure masses. It demonstrates the efficiency of the approximation 
to track possible signatures of substructures. 

We have evaluated these two fields and studied their properties for different lens 
configurations modelled either through massive dark matter halos from a cosmological 
N-body simulation, or via toy models of Monte Carlo distribution of substructures 
embedded in a triaxial Hernquist potential. 

As expected, the angular power spectra of these two fields tend to have larger 
values for larger harmonic numbers when substructures are accounted for and they 
can be approximated by power-laws, whose values are fitted as a function of the profile 
and the distribution of the substructures. 
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1 INTRODUCTION 

The cold dark matter (CDM) paradigm (Cole et al. 2005 
and references therein) has led to a successful explanation of 
the large-scale structure in the galaxy distribution on scales 
0.02 < k < 0.15h Mpc -1 . The CDM power spectrum on 
these scales derived from large redshift surveys such as, for 
instance, the Anglo- Australian 2-degree Field Galaxy Red- 
shift Survey (2dFGRS), is also consistent with the Lyman-Q 
forest data in the redshift range 2 < z < 4 (Croft et al. 2002; 
Viel et al. 2003; Viel, Haehnelt & Springel 2004). 

In spite of these impressive successes, there are still dis- 
crepancies between simulations and observations on scales 
< 1 Mpc, extensively discussed in the recent literature. We 
may mention the sharp central density cusp predicted by 
simulations in dark matter halos and confirmed by the rota- 
tion curves of low surface brightness galaxies (de Blok et al. 
2001) or in bright spiral galaxies (Palunas & Williams 2000; 
Salucci & Burkert 2000; Gentile et al. 2004). Moreover, deep 
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surveys (z > 1 — 2), such as the Las Campanas Infrared Sur- 
vey, HST Deep Field North and Gemini Deep Deep Survey 
(GDDS) are revealing an excess of massive early-type galax- 
ies undergoing "top-down" assembly with high inferred spe- 
cific star formation rates relative to predictions of the hier- 
archical scenario (Glazebrook et al. 2004; Cimatti, Daddi & 
Renzini 2006). 

One problem that requires closer examination concerns 
the large number of sub-L* subhalos present in simulations 
but not observed (Kauffmann, White & Guiderdoni 1993; 
Moore et al. 1999; Klypin et al. 1999). This is the case of 
our Galaxy or M31, although there is mounting evidence for 
a large number of very low mass dwarfs (Belokurov et al. 
2006). However, it is still unclear whether the CDM model 
needs to be modified to include self-interacting (Spergel 
& Steinhardt 2000) or warm dark matter (Bode, Ostriker 
& Turok 2001; Colm, Avila-Reese & Valenzuela 2000) or 
whether new physical mechanisms can dispel such discrep- 
ancies with the observations. For instance, gas cooling can 
be partly prevented by photoionization process which may 



2 S. Peirani et al. 



inhibit star formation in the majority of subhalos (Bullock, 
Kravtsov & Weinberg 2001). 

This "missing satellite problem" remains an ideal 
framework to test cosmological models. During the past 
years, different methods have been employed in order to 
study the gravitational potential of groups or clusters of 
galaxies, for instance through their X-ray lines emission of 
hot gas in the intra-cluster medium or through lensing con- 
siderations. However, while lensing directly probes the mass 
distribution in those objects, the other methods rely more 
often than not on strong hypotheses on the dynamical state 
of the gas and interactions between baryons and dark mat- 
ter. For example, the gas is supposed to be in hydrostatical 
equilibrium in the gravitational potential well created by 
dark matter halo, while spherical symmetry is assumed. In 
this paper, we study the effects induced by substructures on 
the deflection potential of dark matter halos in the strong 
lensing regime. The presence of substructures follows from 
the capture of small satellites which have not yet been dis- 
rupted by tidal forces and/or suggests that the relaxation of 
halos is not totally finished. 

Ray-tracing through N-body cosmological simulations 
suggest that substructures should have a significant impact 
on the formation of giant arcs. At clusters of galaxies scales, 
some results indicate that lensing optical depths can be en- 
hanced (Bartelmann, Steinmetz & Weiss 1995; Fedeli et al. 
2006; Horesh et al. 2005; Meneghetti et al. 2007a) whereas 
some other recent studies suggest that the impact on arc oc- 
curence frequency should only be mild (Hennawi et al. 2007) . 
On the other hand, the presence of substructures changes the 
properties of strongly lensed images to a point that could 
lead to misleading inferred cluster mass properties if not 
properly accounted for (Meneghetti et al. 2007b). Assum- 
ing a one-to-one association of cluster subhalos in the mass 
range 10 11 — 1O 12 ' 5 M0 and galaxies, Natarajan, De Lucia & 
Springel (2007) used weak lensing techniques and found the 
fraction of mass in such subhalos to account for 10 — 20% of 
the total cluster mass, thus in good agreement with predic- 
tions from simulations (Moore et al. 1999). 

Likewise, at the scales of galaxies, strong lensing events 
have long been involving multiple quasars for which the im- 
possibility of resolving the size or shape of the lensed im- 
ages brought the attention toward flux ratios of conjugate 
images as a probe of subtructures. Depart from flux ratios 
expectations from a simple elliptically symmetric potential 
is often interpreted as a signature for local potential per- 
turbations by substructures (Bradac et al. 2002; Dalai & 
Kochanek 2002; Bradac et al. 2004; Kochanek & Dalai 2004; 
Amara et al. 2006). It is unclear whether anomalous flux ra- 
tios actually probe "missing satellites" (Keeton, Gaudi & 
Petters 2003; Mao et al. 2004; Maccio et al. 2006). Due to 
the small source size in the case of lensed QSOs, sensitivity 
to microlensing events due to stars in the lens galaxy makes 
the interpretation less obvious. Astrometric perturbations 
of multiple quasars have also been considered (Chen et al. 
2007) despite substantial observational limitations. 

Presumably the best way out would be to consider ex- 
tended sources like QSOs observed in VLBI or lensed galax- 
ies that will be sensitive to a narrower range of scales for the 
pertubing potential and thus easier to interprete. New meth- 
ods for inverting potential corrections that needed on top of 
a smooth distribution were proposed (Koopmans 2005; Suyu 



& Blandford 2006) but are not guaranteed to converge in all 
practical cases and seem to depend on the starting smooth 
distribution. 

One interesting alternative approach is to treat all devi- 
ations from a circularly symmetrical potential as small per- 
turbations (Alard 2007, 2008) defining the location where 
multiple extended images will form. Two perturbative fields, 
fi(9) and dfo(0)/d8, can then be defined to characterize de- 
flection potential of lenses as a function of the azimuthal 
angle 6, near the Einstein radius. They respectively repre- 
sent the radial and azimuthal derivative of the perturbated 
potential (see Eq. ([10} below). Alard (2008) showed that 
these two fields have specific properties when one substruc- 
ture of mass ~ 1% of the total mass is positioned near the 
critical lines. For instance, the ratio of their angular power 
spectra at harmonic number n is nearly 1. We will investi- 
gate the detailed properties of these perturbative fields by 
considering more realistic lenses such as dark matter halos 
extracted from cosmological simulations. In order to control 
all the free parameters (mass fraction, and shapes of sub- 
tructures for instance) and to study their relative impact on 
arc formation, we will also generate different families of toy 
halos. 

This paper is organized as follows: in section 2 we 
present our lensing modelling; section 3 first sketches the 
pertubative lens solution and applies it to our simulated 
lenses for validation against a ray tracing algorithm; section 
4 presents our main results on the statistics of perturbations, 
while the last section wraps up. 



2 NUMERICAL MODELLING 
2.1 Lens model 

Halos formed in cosmological simulations tend to be cen- 
trally cuspy (p ~ r_1 ) an d are generally not spherical, but 
have an triaxial shape. The triaxiality of these potential 
lenses is expected to increase significantly the number of 
arcs relative to spherical models (Oguri, Lee & Suto 2003 
and references therein), and must be taken into account in 
numerical models. Thus, apart from dark matter halos ex- 
tracted from cosmological simulations, we consider in this 
work typical lenses modelled by of a dark matter halo of to- 
tal mass M = 1O 14 M0 with a generalized Hernquist density 
profile (Hernquist 1990): 
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where R s is the value of the scale radius, R a triaxial radius 
defined by 
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and c/a and b/a the minonmajor and intermediate:major 
axis ratio respectively. We decided to use an Hernquist pro- 
file for practical reasons. However, for direct comparison 
with common descriptions of halos from cosmological simu- 
lation in the literature, the Hernquist profile is related to an 
NFW profile (Navarro, Frenk & White, 1996; 1997) with the 
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same dark matter mass within the virial radius r^ocEl More- 
over, we also impose that the two profiles are identical in the 
inner part (< R 3 ) which can be achieved by using relation 
(2) between R 3 and the NFW scale radius r 3 in Springel, 
Di Matteo & Hernquist (2005). By convention, we use the 
concentration parameter Chost = T 200 /r s in the following to 
characterize the density profile of our lenses. For example, 
for typical lens at a redshift 2 = 0.2, we use R s = 223 kpc 
which corresponds to a NFW profile with Chost = 8.0 (or 
equivalently r2oo = 957 kpc, r s = f f9 kpc) and is consistent 
with values found in previous cosmological N-body simu- 
lations at the specific redshift and in the framework a the 
ACDM cosmology (Bullock et al. 200f ; Dolag, Bartelmann 
& Perrotta 2004). 

Axis ratios of each lens are randomly determined follow- 
ing Shaw et al. 2006: b/a = 0.81 7±0.098, c/b = 0.867±0.067 
and c/a — 0.707± 0.095. These values are in good agreement 
with previous findings from cosmological simulations (see for 
instance Warren et al. 1992; Cole & Lacey 1996, Kasun & 
Evrad 2005). Finally, it is worth mentioning that each halo 
is made of 15 x 10 6 particles corresponding to a mass reso- 
lution of 6.67 x 1O 6 M0. However, we impose a troncation at 
a radius of value 5 Mpc. 

2.2 Substructures model 

2.2.1 Mass function 

Numerical N-body simulations show that dark matter halos 
contain a large number of self-bound substructures, which 
correspond to about 10-20% of their total mass (Moore et 
al. 1999). In the following, the number of substructures N au b 
in the mass range m - m+dm is assumed to obey (Moore 
et al. 1999; Stoehr et al. 2003) 

dA/sub = — r^-dm. (3) 

The normalization constant A is calculated by requiring the 
total mass in the clumps to be 15% of the halo mass and by 
assuming subhalos masses in the range 10 s - 5 x 10 12 Mq. 
The minimum number of particles in the substructures is 
about 15, while the more massive ones have 45, 000 particles. 

2.2.2 Radial distribution 

Substructures are distributed according to the (normalized) 
probability distribution p(r)d 3 r = (p(r)/Mh)d 3 r, where p(r) 
is assumed to have an Hernquist profile of concentration 
C su b, which yields the probability to find a clump at a dis- 
tance r within the volume element d 3 r. While the abundance 
of subtructures in halos of different masses has recently been 
extensively quantified in cosmological simulation (see for in- 
stance Vale & Ostriker 2004; Kravtsov et al. 2004; van den 
Bosh et al. 2007), their radial distribution is less understood. 
However, some studies seem to suggest their radial distribu- 
tion is significantly less concentrated than that of the host 
halo (Ghigna et al. 1998, 2000; Colm et al. 1999; Springel 
et al. 2001; De Lucia et al. 2004; Gao et al. 2004; Nagai 
& Kravtsov 2005; Maccio et al. 2006). We will use either 

1 ^200 defines the sphere within which the mean density is equal 
to 200 times the critical density. 



Csub = 5.0 in good agreement with those past investiga- 
tions, or Csub = Chost for comparison. 

2.2.3 Density profiles and alignment 

The stripping process caused by tidal forces seems to reduce 
the density of a clump at all radii and, in particular, in the 
central regions, producing a density profile with a central 
core (Hayashi et al. 2003). This process was further con- 
firmed by simulations which found that the inner structure 
of subhalos are better described by density profiles shallower 
than NFW (Stoehr et al. 2003). However, other simulations 
seem to indicate that the central regions of clumps are well- 
represented by power law density profiles, which remain un- 
modified even after important tidal stripping (Kazantzidis 
et al. 2004a). This effect may be enhanced when star for- 
mation is taken into account since dissipation of the gas 
(from radiative cooling process) and subsequent star forma- 
tion lead to a steeper dark matter density profile due to 
adiabatic contraction. To test the importance of these dif- 
ferences from the point of view of arcs formation, we allow 
for both of these possibilities: we simulate halos with clumps 
having a central core p(r) oc l/(ro + r) 2 , where ro defines a 
core radius, or a central cusp (Hernquist profile) and study 
how our resulting arcs would change from one option to the 
other. It is worth mentioning that each subhalo concentra- 
tion parameter is obtained using relation (13) in Dolag et 
al. (2004) within the ACDM cosmology. However, to avoid 
spurious effects due to the lack of resolution, all subhalos 
represented by less than 200 particles will have a concentra- 
tion parameter value corresponding to that of an halo made 
of exactly 200 particles (i.e m = 1.33 x 10 9 M Q ). For core 
profiles, we follow Hayashi et al. (2003) and take a core of 
size ro ~ r s . 

Finally, recent cosmological simulations suggest that 
subhalos tend to be more spherical than their host (Pereira 
et al. 2008; Knebe et al. 2008) and this effect can also be 
enhanced if halos are formed in simulations with gas cool- 
ing (Kazantzidis et al. 2004b). Moreover, the distribution 
of the major axes of substructures seems to be anisotropic, 
the majority of which pointing towards the center of mass 
of the host (Aubert, Pichon & Colombi 2004; Pereira et al. 
2008). Although shapes and orientations of subhalos provide 
important constraints on structure formation and evolution, 
we have not studied their relative influence in this work. We 
reasonably think that modelling substructures by either tri- 
axial shapes of spherical shapes won't lead to any significant 
differences in our results. 



2.3 Lens samples 

Table 1 summarizes our different samples of lenses. From a 
statistical point of view, each sample involves one hundred 
realizations of halos following the above described method- 
ology. Lenses are assumed to be at a typical redshift z — 0.2. 
They share a common total mass of 10 14 Mq and are thus 
described by an Hernquist profile of concentration C = 8. 
The sample A represents our reference catalogue in which 
all halos have no substructure. For each of them, we intro- 
duce a fraction of subtructures (F su \,) by removing some 
background particles so that both the total mass and the 
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Figure 1. Elliptical luminosity contour of a source with r]o = 
0.2. The three dashed lines correspond to isophotes defined by 
I(Rl) = 0.6 I{Ri) = 0.2 and I(R 3 ) = 0.01 



density profile of the initial halo are conserved. These halos 
are classified in catalogues B and C according the definition 
of inner density profile (IP) of clumps. For example, each 
halo from samples Bl and B2 have 15% of substructures 
with an inner profile represented by a cusp (Hernquist pro- 
file). Their radial distribution (RD) within the halo is left 
as a free parameter. We consider two possibilities, C su b = 5 
as suggested by numerical simulations, and C 3U b ~ 8, which 
is the concentration parameter of the whole halo. Finally, 
lenses catalogues CI and C2 have substructures represented 
by a core profile with C 3U b = 5 and C SU b = 8 respectively. 

In addition, lensing efficiency depends on the relative 
distance between lenses and sources. The efficiency of the 
lens is scaled by the critical density 
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where D s , D s and Da s are the angular diameter distances 
between the observer and the source, between the observer 
and deflecting lens and between the deflector and the source 
respectively. When the surface mass density in the lens ex- 
ceeds the critical value, multiple imaging occurs. In order 
to account this effect which, for a given lensing halo, im- 
plies a different Einstein radius for a different source red- 
shift, we consider the redshift distribution of sources taken 
from the COSMOS sample of faint galaxies detected in the 
ACS/F814W band (Leauthaud et al. 2007). It is well repre- 
sented by the following expression 
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with z = 0.345 and a = 3.89 (Gavazzi et al. 2007). 
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3 NUMERICAL VALIDATION OF THE 
PERTURB ATIVE SOLUTION 

In this section, we take advantage of the large sample of 
mock lenses described in § [2] to assess the validity of the 
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Csub — 8 



Table 1. Samples of lenses (see details in the text). 

perturbative method developed in (Alard 2007). After a brief 
presentation of the basic idea in ol l3.ll we compare the ability 
of this simplified procedure to reproduce multiple images 
lensed by complex potentials as compared to a direct ray- 
tracing method (§ I3.3f) and define the validity range of the 
perturbative approach. 

3.1 The perturbative approach 

For the sake of coherence, let us sketch the motivation 
behind the perturbative lens method (Alard 2007) used 
througout this paper. The general lens equation, relating the 
position of an image on the lens plane to that of the source 
on the source plane can be written in polar coordinates as 



/ d(j>\ fldcj>\ 



(6) 



0, 



(7) 



where r s is the source position, and r, u,- and ue are the 
radial distance, radial direction and orthoradial direction re- 
spectively. Here </>(r, 8) is the projected potential. Let us now 
consider a lens with a projected density, S(r), presenting cir- 
cular symmetry, centered at the origin, and dense enough to 
reach critical density at the Einstein radius, Re - Under these 
assumptions, the image by the lens of a point source placed 
at the origin is a perfect ring, and equation (]§} becomes: 

r _ 

dr 

where the potential, <f>o, is a function of r only, and the 
zero subscript refers to the unperturbed solution. The basics 
ideas of the perturbative approach is to expand equation Q 
by introducing i) small displacements of the source from the 
origin and ii) non-circular perturbation of the potential, ?/> 
which can be described by: 

r s = er s , and </> = 4> + , (8) 

where e is small number: e <C 1. To obtain image positions 
(r,8) by solving equation ([6} directly, may prove to be analy- 
ically impossible in the general case. It is then easier to find 
perturbative solution by inserting equation (JSJ into equation 
((6|. For convenience, we re-scale the coordinate system so 
that the Einstein radius is equal to unity. The response to 
the perturbation on r may then be written as 

r = 1 + edr , (9) 

which defines dr(9), the azimuthally dependant enveloppe 
of the relative deflection. Using Equation (JHJl, the Taylor 
expansion of <f> is 



= 4>o + e4> = J2[C n + e/„(0)](r - 1)' 



(10) 
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where: 



, and 

n! L dr" J r=i 



(ii) 



(12) 



Finally, inserting equations and (JTOJ into equation (JSJ) 
leads to: 



la t \ dfo 

r s = (K2flr — ,/i)u r jf U9 • 

dp 



(13) 



with K2 = 1 — 2C2. This equation corresponds to Eq. (8) in 
Alard (2007). 



3.2 Morphological effects VS astrometric 
distortions 

As demonstrated in Alard (2008), the morphology of arcs is 
very sensitive to small perturbators such as substructures in 
the main halo. The effect of substructures on morphologi- 
cal features like for instance the size of an image is typically 
much larger than pure astrometric distortions. Indeed, astro- 
metric distortions are only of the order of the substructure 
field, which in many case is smaller than the PSF size of 
the instrument used, and represents then an observational 
challenge to measure. Note that a perturbative theory of as- 
trometric distorsions was already considered by Kochanek 
et al. (2001) and Yoo et al. (2005; 2006). Note also that the 
former perturbative approach is limited to astrometric ef- 
fects, the relevant theory does not describe image formation, 
and thus cannot predict morphological effects. When image 
morphology is considered, effects are an order of magnitude 
larger than astrometric effects. For instance, let us consider 
the examples displayed in Fig.l and Fig. 2 of Alard (2008). 
Whereas a giant arc is obtained from an unperturbed ellip- 
tical lens cusp caustic (Fig. 1), the introduction near the 
Einstein radius of a substructure of only 1 % of the main 
halo mass breaks the arcs into 3 sub-images (Fig. 2). The 
detection of such effects should not require a particularly 
good resolution, since the amplitude of the effect is a frac- 
tion of the arc size, which is typically several time the PSF 
size. 

Specifically, let us consider the following example: we 
take a sub-critical configuration for an elliptical lens and 
evaluate the modification of the image size due to the per- 
turbation by a substructure field. Let us assume that the 
ellipticity of the lens is aligned with the axis system and 
both source and the substructure are placed on the X-axis. 
In this configuration, the size of the central image will be 
perturbed by the substructure field, and this will be the 
observable effect. The main interest of this example is the 
simplicity of the calculations (linear local description of the 
field) and the simple description of the effect (reduction of 
image size). Let us now define some useful quantities: the 
local slope of the unperturbed field ao, the angular size of 
the unperturbed image Aft) and the source size Ro. The sub- 
structure parameters are its mass mp (in unit of the main 
halo mass) and its distance dr to the Einstein ring. For cir- 
cular sources, the local slope, image size and source radius 



are related by the following relation in the unperturbed case: 
A6 ° " ' (14) 
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The slope perturbation introduced by the substructure da 
is (See Alard 2008, Sec. 2): 

mp 



da ■ 



dr 



(15) 



The substructure field modifies the images size according to: 
(a + da) (AOo + d9) = R , (16) 

which gives: 

da 2 mp Ro 



= Aft 



0- 



dr Qq 



(17) 



For an elliptical isothermal potential with ellipticity param- 
eter 77, the field d/o(0)/d0 reads (see Alard 2007, 2008): 

dfo 



dO 



= r\ sin 29 — T) sin ( 



(18) 



Thus, ao which is the field derivative in is: ao — rj , leading 
to: 



dr rj z 



(19) 



The former equation evaluates the angular perturbation of 
the image size by the substructure. To obtain the corre- 
sponding perturbation on the image length dS, we have to 
multiply by the Einstein Radius: 



dS~R E \da\ =2R E ^^ 
dr n-* 



(20) 



Note that astrometric effects are identical to the effects of 
the substructure on the /1 field for circular sources. Thus, 
the astrometric effect dA is of the order (See Alard 2008, 
Sec. 2): 

nip 



dA ~ Re- 



dr 



(21) 



Consequently, the ratio between morphological and astro- 
metric effects here is: 

dS 2R 



dA 



T 



(22) 



For arcs, the source size (source diameter= 2Ro), and the 
parameter rj have the same typical scale, which gives: n = 
0.1, and Rq = 0.05; then: 



dS_ 

dA 



10. 



(23) 



This means that morphological effects are 10 times larger 
than astrometric effects. This point is critical, since it really 
makes the effect of substructure observable. To illustrate 
this latter point, let's consider some numerical values for 
typical galaxies. For a Milky way like galaxy, we have a 
mass M ~ 6 10 Mq. The Einstein radius in arcsec is given 
by: 



Re 



M 



10 12 M, 



= 1.4 arcsec , 







(24) 



which gives the size of the image perturbation, dS 
Re \da\: 



dS ~ 2.8^7^- -^7- arcsec . 
dr rj 2 



(25) 
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For a perturbator with 0.5 % mass of the galaxy, taking 
typical scales, dr = 0.1, Ro = 0.1, and r\ — 0.1, we obtain: 



dS ~ 1.4 arcsec . 



(26) 



Such effects should be within the reach of spatial instru- 
ments such as DUNE or SNAP. On the contrary, astrometric 
effects are typically 10 times smaller dA ~ 0.14 arcsec and 
should be then much more difficult to detect. Consequently 
morphological effects are definitely our best hopes to detect 
substructures. 

Clearly, other configurations will also allow the mea- 
surement of morphological effect of substructures, though 
in general the calculations will be a little bit more compli- 
cated, but the effects will more or less be of the same order, 
since they are related to the amplitude of the perturbation 
of the dfo(8)/d9 field by the substructure. 

3.3 Reconstruction of images 

One interesting feature of the perturbative method is to pro- 
vide a framework for the reconstructions of images. By first 
defining an elliptical source centered on position (xo,J/o), 
with a characteristic size Ro, ellipticity y/2r/o, and inclina- 
tion of the main axis #o such that, 

Ro = {1 - riQCosOo)(x a - x ) 2 

+ 2rj sin26»o(a; s - x Q )[y s - y ) 



situations. In particular, we use lenses either from cosmolog- 
ical simulations or from toy models presented in section [2] 

For direct comparison between arc reconstructions pre- 
dicted by the perturbative method and theoritical ones, we 
use a ray-tracing method. Part of our investigations indeed 
makes use of the Smooth Particle Lensing technique (SPL), 
described in details in Aubert, Amara & Metcalf (2007) 
and summarized in this section. SPL has been developed 
to compute the gravitational lensing signal produced by an 
arbitrary distribution of particles, such as the ones provided 
by numerical simulations. It describes particles as individ- 
ual light deflectors where their surface density is arbitrarily 
chosen to be 2D Gaussian. This choice makes it possible to 
compute the analytical corresponding deflection potential, 
given by: 



+ (1 + ??ocos6>o)(2/ s - y ) 



(27) 



one can express the equations of the image contours using 
Eq. (J33J: (see Alard (2007) for more details) 



f 

K 2 



r | yjRjS-(l-^)(df /der 



S = l-r ?0 cos26», 



(28) 



This equation corresponds to Eq. (15) in Alard (2007). The 
functional fi is defined to take into account the effect of the 
translation of the source by the vector ro = (xo,yo)- 



fi = fi + xo cos + y Q sin ( 



for 



0,1. 



(29) 



As emphazised in Alard (2007), the image contours are only 
governed by the two fields, fi(0) and dfo(0)/d8 , which con- 
tain all the information on the deflection potential at this 
order in the perturbation. For instance, the two first terms 
in the bracket of Eq. (I28|l give informations on the mean po- 
sition of the two contour lines while the last term provides 
informations on the image's width along the radial direction 
as well as a condition for image formation. Therefore, the 
characterization of these two fields represents a simple and 
efficient way to track possible signatures of the deflection po- 
tential induced by substructures, so long as the perturbative 
framework holds, as we will illustrate below. Alard (2007) 
already implemented the method with a lens described by a 
NFW profile yielding an analytical solution for the projected 
potential profile. In this section, we illustrate and validate 
the method while considering more complicated and realistic 



W = ^(log(^)-2Ei(-^)), 



47rE c v '~ M C r 4y "~' v 2a 2 "' (30) 

where Ei(x) = — f°° exp(— x)/xdx, m p is the mass of the 
particle, a its extent and E c the critical density. From the 
deflection potential, expressions for the deflection angles a, 
the shear components 7 and the convergence k can be easily 
recovered (see Aubert, Amara & Metcalf (2007) for more 
details). Knowing the lensing properties of a single particle, 
one can recover the full signal at a given ray's position on 
the sky by summing the contributions of all the individual 
deflectors: 



i i 

K ( r ) = J2 Ki ^> 7(r) = ^Z^W' 



(31) 
(32) 



where, e.g. %(r) is the contribution of the i-th particle to 
the shear at ray's position r. This summations are per- 
formed efficiently by means of 2D-Tree based algorithm, in 
the spirit of N-body calculations. The tree calculations are 
restricted to monopolar approximations where an opening 
angle of 0.5 — 0.7 is found to give results accurate at the 
percent level on analytical models. Finally, Aubert, Amara 
& Metcalf (2007) found that an adaptative resolution (i.e. 
an adaptative extent a for the particles) provides a signifi- 
cant improvement in the calculations in terms of accuracy. 
For this reason, the smoothing a depends on the rays lo- 
cation: particles shrink in high density regions in order to 
increase the resolution while they expand in low-density re- 
gions, smoothing the signal in undersampled areas. For all 
simulations, we use 1024 x 1024 rays within a square of size 
2 x 2Re, an opening angle of 0.7 and N a = 256 (where 
N a is the number of particles over which the smoothing is 
applied) . 

3.3.1 Lenses from the toy model 

We present in this section characteristic examples of arc re- 
construction. Three lenses Lo, Li and L2 belonging to sam- 
ples A, B2 and C2 respectively are considered. They have a 
common mass, density profile, axis ratios and random orien- 
tation in 3D space. They only differ via the presence or not 
of substructures as well as via the inner density profile of 
substructures: Lo has no substructure whereas Li has sub- 
structures with a central cusp while a core describe the inner 
density profiles of substructures in L2. In the present case, 
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Einstein radius unit 

Figure 2. Projected density maps of lenses modelled by toy halos and their associated arcs reconstructions when considering an elliptical 
source contour with a Gaussian luminosity. Panels (a) and (e) show respectively the projected density map of the lens Lq and the resulting 
image obtained from ray-tracing when the source is placed at the origin. Both red and light blue dashed lines represent the Einstein 
radius. The solid white line is the predicted arc reconstruction for isophotes O.OlTmax. Similar plots are shown in panels (b) and (f) but 
for L3. Columns 3 and 4 show arc reconstructions for the lens Li. The solution from ray-tracing is given in the panel (c) while the other 
last panels present a direct comparison between the isophotal contour (red) and the arc reconstruction (blue) for O.Olimax (g), 0.2/ m ax 
(d) and 0.6/ m ax (h). 



the source is at a redshift z ~ 2.9, has an elliptical contour 
with 770 = 0.2 and a radius Ro ~ 0.05-Re, characterizing a 
Gaussian luminosity profile (see figure [l}. 

For the arc reconstructions presented below, we shall 
consider 3 different radii Ri, R2 and R3 corresponding to 
3 specific isophotal contours defined by I(Ri) = 0.6 7 max , 
I{Ri) = 0.2 J ma x and I(R 3 ) = 0.01 / max (see figured}. 

In Figure @, we show the projected mass density of 
lenses Lo and L\ near the Einstein radius, the image's so- 
lution obtained from ray-tracing and the contours predicted 
by the perturbative method when the source is placed at 
the origin. When no substructure is considered, both pro- 
jected density and potential are nearly elliptical. As ex- 
pected, we obtain four distinct arcs in a cross configuration. 
The predicted arcs reconstruction are in good agreement 
with the numerical solution obtained via the ray-tracing al- 
gorithm. From a theoretical point of view, it is easy to show 
that the functions fi(0) and dfo(9)/d6 are proportional to 
oc cos(2(9 + ip) and oc sin(2# + ip) respectively. These func- 
tionnal form are recovered in our experiment and shown in 
Figure ©. 

When substructures are present, the shape of images is 
significantly altered. First, we notice that the positions of 
substructure tend to break the ellipticity of the halo center. 
Thus, it is not surprising that the shape of the image is ap- 
proaching a ring in that case. Moreover, it is interesting to 
see that the position of one substructure (at the top left of 
the figure) is exactly at the Einstein radius. This produces 
an alteration of the luminosity, while the effects is more vi- 



olent when substructures present a cups profile. This effect 
can be clearly seen when comparing the perturbative fields 
fi(9) and df (9)/d9 relative to the lens L\ in Figure ©. 
For instance, we can see two clear bumps in the evolution of 
dfo(9)/d6 . The second one (9 > 2ir/3) is produced by sub- 
structures in the lower right part and induced an alteration 
of the luminosity again. 

To estimate the systematic error between the theoreti- 
cal contours provided by the ray-tracing and those predicted 
from equation l|28p , we use a simple procedure with a low 
computational cost. First, each predicted contour is divided 
into a sample of N points. Each of them is defined by po- 
lar coordinates (ri, 9%) which coincide with a luminosity 
value of the image corresponding to an unique radius -Ri 
in the source frame. By using relation (|28[) . we then com- 
pute 1 + dr(Ri) which gives the image contour radius of 
the isophotal contour I(Ri) of the source in Einstein radius 
unit. By defining 1 + dri the radial distance of point i, the 
mean error err (in Einstein radius unit) is then computed 
byerr = £f \dr(Ri) - dn\/N. 

For illustration, we have estimated the mean error 
reached for the lenses L\ and Li using three luminosity 
contours source (see fig. [l}. For L\, the mean errors are 
respectively 0.67%, 0.71%, 0.95% of the Einstein radius for 
isophotes 0.6, 0.2 and 0.017 ma x respectively, while we obtain 
0.74%, 0.86% and 1.04% R E for lens L 2 and for same lumi- 
nosities. We have also studied how the mean error evolves 
for random positions of the source inside an area limited by 
the caustic lines. To do that, we have used the lens L\ and 
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Figure 3. Variations of the fields fx and d/o(£>)/d£> as a function 
of 8 for the lens Lq (panels a and b), lens L\ (panels c and d) 
and lens L2 (panels e and f). 



have studied 100 realizations with different impact parame- 
ters. We found err = (1.01 ± 0.12)%_Re for isophotes equal 
to 0.01/ max . 
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Figure 4. Projected density map (panel a) of a lens modelled 
by a dark matter halo extracted from the N-body simulation and 
the associated arc reconstructions. The solution from the ray- 
tracing is plotted in the panel (b) with the Einstein radius (blue 
line) and the predicted arc reconstruction for isophote O.OlTmax- 
The second line compares the isophote contours 0.01/max (panel 
c) and 0.27 max (panel d) represented by the blue lines with the 
predicted contour (red lines). The same results are presented in 
the panel e for isophote 0.6/ max . Finally, variations of fi(8) (red 
line) and d/o(9)/d6> (black line) are plotted in the panel (f)). 



3.3.2 Lenses from cosmological simulation halos 

In this section, lenses are modelled by dark matter ha- 
los extracted from a cosmological simulation of the Pro jet 
HORIZOI0. The simulation was run with Gadget-2 (Springel 
2005) for a ACDM universe with Q M = 0.3, fi A = 0.7, 
Q B = 0.045, Ho = 70 km/s/Mpc, a s = 0.92 in a periodic 
box of 20 /i _1 Mpc. We use 512 3 particles corresponding to 
a mass resolution of m par t. — 4 x 1O 6 M0 and a spatial reso- 
lution of 2 kpc (physical). Initial conditions has been gener- 
ated from the MPgrafic code (Prunet et al. 2008), a parallel 
(MPI) version of Graf ic (Bertschinger 2001). In this simu- 
lation, we selected two regions. In the first one, the lens is a 
typical halo of total mass 2.6 x 1O 13 M at redshift 2 = 0.5. 
The source is at 2 = 1.2, assumed to be elliptical (770 = 0.2) 

2 http://www.projet-horizon.fr/ 



with Rq — 0.05-Re and placed near a caustic in order to ob- 
tain a giant arc. Fig. Q shows the projected density of the 
lens near the Einstein radius; both the ray-tracing solution 
and the predicted contours by the perturbative method are 
shown. Here again, the three different contours are well re- 
constructed since the error are 0.76%, 0.83% and 0.91% Re 
for isophotes 0.6, 0.2 and 0.01/ max respectively. For illustra- 
tion, the angular variation of the perturbative fields are also 
representated in the figure @. 

The second example is a lens modelled by another 
halo from the same N-body simulation. Its total mass is 
6.6 x 1O 13 M0 at 2 = 0.5. This is an extreme case since a sig- 
nificant number of substructures are still falling toward the 
center of the host halo which suggests that the dynamical 
relaxation is still operating. This strongly affects the po- 
tential and the perturbative fields (see figure [5}. However, 
mean errors remain small of the order of 1.11%, 1.20% and 
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Figure 5. same as Figure f4] for another dark matter halo ex- 
tracted from the N-body simulation. 

1.26% Re for isophotes 0.6, 0.2 and 0.01I ma x respectively, 
which proves the accuracy of the method to deal with more 
complex systems. We may reasonably think that lenses in 
our different samples at z = 0.2 tend to be more relaxed 
that the present configuration, and, consequently, the mean 
errors in our arc reconstruction should be less pronounced. 



each couple of parameters (r/o, Go)- The proper solution cor- 
responding to the true parameter (r)o,9o) has minimal prop- 
erties. Consider for instance a circular solution (770 = 0), if 
the inversion formula is used with additional Fourier 

terms with order n > 2 will appear in the inversion formulae. 
Thus, it is clear that minimizing the power in higher order 
Fourier modes is a criteria that will allow to select the best 
solution when exploring the plane (770, #o)- This criteria has 
also a very interesting property, considering that power at 
order n > 2 usually reveal the presence of substructures (see 
Table 3 for instance), the solution with minimum power at 
higher order is also the one that puts the more robust con- 
straint on the presence of substructure. Thus the elliptical 
inversion can be performed by exploring the plane (770,60) 
in a given parameter range, computing the corresponding 
fields fi(0) and dfo(9)/d9 and their Fourier expansion, and 
selecting the solution with minimum power at n > 2. For non 
elliptical sources, one can use the general inversion method 
presented in Alard (2008). This inversion method remaps 
the images to the source plane using local fields models (ba- 
sically the scale of the images). The solution is selected by 
requiring maximum similaritiy of the images in the source 
plane. Image similarity is evaluated by comparing the im- 
age moments up to order N. Provided the number of image 
moments equations exceed the number of model parame- 
ters, the system is closed and has a definite solution. Note 
that the local models may be replaced with general Fourier 
expansion in the interval < 9 < n, but in this case, the 
additional constraint that no image are formed in dark areas 
must be implemented (Diego et al. 2005). In the perturba- 
tive approach this requirement can be reduced to dfo(9)/d9 
> Rc in dark areas, where Rc is the radius of the smallest 
circular contour that contains the source. We may therefore 
assume for now that observational data may be inverted, 
and that an observationnal survey of arcs should provide 
us with a statistical distribution of the perturbatives fields. 
Hence we may use our different samples of halos presented 
in paragraph 12.31 in order to measure the relative influence 
on arc formation of the different free parameters such as 
the inner profile of substructures or their radial distribution 
within the host halo. To conduct this general analysis the 
fields will be represented by Fourier models, due to the direct 
correspondance between Fourier models of the fields and the 
multipolar expansion of the potential at r = 1 (Alard 2008). 



4 FOURIER SERIES EXPANSION 
4.1 Motivation 



Given Eq. ((28}, it is straightforward, for a given lensed im- 
age and an assumed underlying spherical lense and elliptical 
source, to invert it for f\(9) and dfo(9)/d9 as 



±S 



d t " tL=Jml/S - 4[dr+(6) - dr~(9)f (33) 

1 - T7n V 



where S, 9 and K2 are given by Eq. (|28|) . This inversion for- 
mula depends explicitly on the source parameters, (770,6*0) 
which are unknown. However using Eqs ( 1331 I34p it is pos- 
sible to compute the two functions fi(0) and dfo(9)/d9 for 



4.2 Results 

The angular functions fi(6) and dfo(9)/d8 can be charac- 
terized by their Fourier expansion: 



dfo(6) 
d9 



ft(9) = j2< a ^ cos ( nd + ^) 



Pi(n) 



((al) }, where i = 0, 1, 



(35) 
(36) 
(37) 



where Pi(n), i = 1,2 correspond to associated power spec- 
tra. We have derived the multipole expansion of fi(9) and 
dfo(9)/d9 for each halo of the different catalogues and we 
focus in the following on the mean amplitudes (a°) and (a* } 
obtained. 
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Table 3. Power spectra of d/o(0)/dfl shown in the second column 
of Figure [3] . 



Tables ([2]) and ([3| respectively summarise the seven first 
orders of the power spectrum of fi(0) and dfo(9)/d9 for the 
3 lenses Lo, L\ and L%. 

When substructures are absent, both harmonic power 
spectra of fi(0) and dfo(9)/d9 are dominated by the second 
order mode, which is characteristic of a projected elliptical 
potential. The situation is totally different when substruc- 
tures are taken into account. First, we notice that first mode 
(n — 1) increase for lenses L\ and L 2 . This is due to the fact 
that we kept the same definition of the mass center between 
the three lenses. The random position of subtructures gener- 
ates a non zero impact parameter which affect the first order 
mode according the relation (129|l . Moreover, since substruc- 
tures tend to break the ellipticity of the halo center in the 
present case, one expects that the second mode decreases. 
However, the most interesting feature is that modes corre- 
sponding to n > 3 increase when substructures are present. 

4.3 Practical limitations 

A fraction of the error is produced by the ray-tracing simu- 
lation as well as the limitation of the considered resolution. 
Let us therefore consider a toy halo with an isothermal pro- 
file: 



Po 

„2 ■ 



(38) 



where po is evaluated so that the mass enclosed inside a ra- 
dius r — 957 kpc is M — 1O 14 M0. An isothermal profile 
is appropriate to estimate systematic error since it leads to 
exact solution with the perturbative method. Here again, we 
have evaluated the quantity err by considering 100 realiza- 
tions with different impact parameters. Sources have circu- 
lar contour with Ro ~ 0.05Re and 15 millions particles have 
been used. For isophotes 0.6I ma x, which is supposed to have 
the higher error values, we have obtained err = (0.30±0.03) 
% Re- Thus, in the following, we will consider that both 
ray-tracing method and the resolution limitation induce to 
a mean error of 0.3 % Re in contours reconstructions. More- 
over, as we shall see in section[4] fi(6) and dfo(9)/d9 can be 



characterized by their multipole expansion and their associ- 
ated power spectrum (see equations [37)) . In Figure (JB), we 
plot the amplitudes (a 9 ) and (a^} as a function of n derived 
from the present experiments. These values for the ampli- 
tude (ct ~ 0.06%-Re) correspond to a noise that we have to 
take into account below. For this reason, we put a confident 
limit to ~ 2. Oct. 

In Aubert et al. 2007, the influence of smoothing 
number of particles N par t and opening angles 9 have been 
extensively investigated on softened isothermal spheres. The 
number of particles and the models used here are similar 
and the parameters used in the current study can be con- 
sidered as the most appropriate considering these previous 
tests (N a = 256, 9 = 0.7, N part = 15 X 10 9 ). For instance 
a re-analysis of the Aubert et al. tests case present an aver- 
age error in the deflection angle of 0.3-0.4 % at the Einstein 
radius. Considering that the profiles and the number of par- 
ticles are similar, the current error estimation is in good 
agreement. 

In the same studies, the (inverse) magnification recon- 
struction was also previously tested and the critical lines 
fluctuates around their theoretical location because of Pois- 
son noise and ray-shooting artefacts. By increasing the num- 
ber of particles and by means of adaptive smoothing, these 
errors can be limited. Again, the re-analysis of these tests 
cases allows an estimation of the error on magnification p, 
of Sfi/fi ~ 0.02/x close to the Einstein's radius. The inverse 
magnification is obtained from the joint calculation of the 
convergence and the shear through the same ray-shooting 
technique. Hence, the error estimation does not rely on some 
propagation procedure but on the effective calculation and 
thereby includes Poisson sampling effects and ray-shoot er- 
rors. 

To finish, it's interesting to determine how the error of 
0.3 % Re in arc reconstruction is reflected in error in the 
image length. Let's consider again the configuration studied 
in Sect. 13.21 and taking equations (|14[) and (|17[) . one obtains: 



\d0\ = 2A0i 



da 
o — 

Qo 



. da Ri 



(i 



(39) 



Taking da = 0.3/100 with the same typical value, Ro — r\ 
0.05, and qq = r\ — 0.1 we have: 



\d9\ = — 
1 100 



(40) 



As in the previous calculation, the actual size of the pertur- 
bation due to the error dE is obtained by direct multiplica- 
tion with the Einstein radius. We take the same numerical 
value, Re = L4 arcsec, thus: 



dE = 0.04 arcsec 



(41) 



Given the amplitude of morphological effects, this noise is 
not a source of concern, but had we considered astrometric 
effects, this noise would become be a real problem. 

4.4 Statistics 

Figures (0 and Q show respectively the variations of the 
mean amplitudes (a°) and (aj,) as a function of n derived 
from our different samples of lenses. As mentioned above, we 
put a confident limit (~ 2. Oct = 0.12%.Re) and we exclude in 
our calculation all amplitudes below. First, when substruc- 
tures are disregarded, power spectra are dominated by the 
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Figure 6. Variations of the mean amplitudes, (a^) (left panel), 
and (a*) (right panel), as a function of the harmonic order, n, 
derived from 100 lenses modelled by an isothermal profile. 



second harmonic which just reflects the fact that our simu- 
lated lenses have a mean ellipticity. However, values of the 
fourth orders appear to be not negligible. This is probably 
due to the 2D projection which can lead to boxy projected 
densities. On the other hand, odd orders are negligible, as 
expected. When substructures are present, we can notice 
that the amplitudes of first harmonic, (oi), have rather high 
values. As emphasized before, this is simply due to the fact 
the position of the lenses center is assumed to be the same 
for all objects. In fact, subtructures tend to modify the po- 
sition of the mass center or, equivalently, tend to generate 
a non zero impact parameter which affect first harmonic co- 
eficient according to the relation (|29[) . However, the most 
interesting and important result is the presence of a tail in 
the power spectra which clearly suggests that the amplitude 
of high order harmonics (n > 3) are not negligible anymore. 
Note that these effects are clearer when the density profile of 
substructures is modelled by a cusp as was the case for sam- 
ple B2. Note that the power spectra of high orders harmonic 
(n > 3) can be fitted by power-laws: 



Po(n > 3) = k n a ° , P 1 (n>3)=k 1 n° 



(42) 



where the relevant fitting parameters are shown in the Ta- 
ble n 

Having estimated the statistical contribution of sub- 
structures to the perturbative fields, and we may now ad- 
dress the problem of defining observational signatures of 
substructures. The observational effects are of two types: 
(i) an effect on the position of images which is controlled by 
the field fi (9) , and (ii) effects on the image morphology, for 
instance its size in the orthoradial direction which depends 
on the structure of the field dfo(9)/d9 . An interesting point 
is that effects of type (i) are directly related to /i (9) , while 
type (ii) effets is not related directly to dfo(6)/d9 but to its 
derivative (at least for images of small extension). Indeed, 
let us consider a source with circular contour; provided the 
image is small enough, a local linear expansion of the field 
dfo(9)/d6 will be sufficient to estimate the image morphol- 
ogy. We make the following field model: 



de 



de 2 



(9 - 0q) . 



(43) 



where 9o is a position angle which should be close to the 
image center. At the edge of the images, we have: 



Table 4. Fit parameter of the statistical distribution of the har- 
monics of the two fields. 



By defining A9 — 9 — 9o, it follows that: 
d 2 /o 



d<9 2 



■A0 = Ro ■ 



(45) 



Thus, (d 2 /o(#)/df? 2 ) -1 scales like the orthoradial image size 
Af? and can be directly related to observational quantities. 
The main difference with effects of types (i), which are re- 
lated to the field rather than the field derivative, is that 
the derivative introduce heavy weights on the higher orders 
of the Fourier serie expansion of the field, where the sub- 
structure contribution is dominant. Indeed, the derivative 
of the Fourier serie introduce a factor n at order n, which 
translates in a factor n 2 on the components of the power 
spectrum. Thus the image morphology (and in particular its 
orthoradial extension for smaller images) will be much more 
sensitive to substructure than the average image position. 
To illustrate this, we have plotted in Figure (|9} the ampli- 
tudes (a°) for lenses samples Bl, B2, CI and C2. Moreover, 
to study the contribution of order n > 3, it is convenient to 
define the following quantities: 

Pi = ((«!) 2 ) + 4x«4) 2 ) 

I> 2 x<(a°) 2 ) 



(46) 



R - dfo 



(44) 



P 2 = 

n>3 

/pT 

Q = VpT 

Table |4] summarize mean values of Q relative to each cat- 
alogue of lenses. As expected, Q has higher values when 
substructure have a cusp profile. The total contribution 
of high order vary between ~ 8 and 18% according the 
model of lenses used, which is quite significant. In clos- 
ing, for the sample C2, mean errors between predicted con- 
tours from the perturbative method and ray-tracing solution 
are err = (0.95 ± 0A7)%R E , err = (0.97 ± 0A8)%R E and 
err = (1.11 ± 0A9)%R B for isophotes 0.6, 0.2 and 0.01/ max 
respectively. 



5 CONCLUSIONS 

The structuration of matter on galactic scales remains a 
privileged framework to test cosmological models. In partic- 
ular, the large number of dark matter subhalos predicted by 
the ACDM cosmology is still a matter of debate. In this pa- 
per, potential signatures of substructures in the strong lens- 
ing regime were considered. This investigation makes use of 
the perturbative solution presented by Alard (2007, 2008), 
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in which small deviations from the "perfect ring" configu- 
ration are treated as perturbations. In this framework, all 
informations on the deflection potential are contained in two 
one-dimensionnal fields, fi(0) and d/o(#)/d# which are re- 
lated to the radial expansion of the perturbed aspherical 
potential. The analysis of the properties of these two fields 
via their harmonic decomposition represents a simple and 
efficient way to track back possible signatures of subtruc- 
tures. For instance, the perturbative method offer a simple 
and clear explanation of the multiplicity of the images. As 
shown in Alard (2008) , the unperturbed image is a long arc 
for a cusp caustic configuration, but when a substructure is 
introduced near the arc, this latter is broken in 3 images. 
This effect is also illustrated in the present work in the ex- 
ample shown in the figure . Basically the breaking of the 
image and the change in image mutiplicity is due to the per- 
turbative terms introduced by the substructure on the field 
d/o(f?) /d(9. In the case of a circular source with radius Ro, 
the image is broken when d/o(#)/df? > _Ro- 

In this paper, lenses were modelled either by dark mat- 
ter halos extracted from cosmological simulations or via toy 
models. The advantage of toy models is to reach a higher 
resolution and to allow us to study the influence of free pa- 
rameters such as the inner profiles of substructures which 
are expected to play a central role here. We have first esti- 
mated the accuracy of the perturbative predictions by com- 
paring the mean error between predicted contours images 
and theoretical ones derived via ray-tracing. We found that 
in general, the relative mean error is ~ 1% of the Einstein ra- 
dius {Re) when different impact parameters configurations 
and different source redshifts are considered. We also found 
that both resolution limitation and ray-tracing procedure 
lead to a systematic error of ~ 0.3%Re- Furthermore, al- 



though the accuracy of this approach for elliptical lenses is 
demonstrated in Alard (2007), we have checked and shown 
in the present study that the method works for the range of 
ellipticities derived from realistic 3D dark matter halos. This 
implies that the numerical evaluation of the coefficient of the 
perturbed potential, fi(9) and d/o(f?)/df? at the Einstein ra- 
dius is accurate enough to carry a statistical investigation. 

We have generated several mock catalogues of lenses 
in which all objects have a total mass of 1O 14 M0 and are 
at redshift z = 0.2. This value is motivated by comparison 
with observational surveys, and is close to where the strong 
lensing efficiency of clusters is the largest for sources z s > 1 
(Li et al. 2005). The first catalogue represents our reference 
sample since all lenses are modelled by dark matter halos 
without substructures. In the other ones, substructures are 
described by either a cusp profile or a core profile. Their ra- 
dial distribution is also a free parameter and we have used 
Csub = Ghost and C su b = 5. Our statistical investigation in- 
volves a Monte-Carlo draw: the ellipticity of host halos, po- 
sition of substructures, sources redshift for instance are ran- 
domly derived according to specific distributions. We found 
that the harmonic power spectra of /i (ff) and d/o (8) /d8 tend 
to develop a tail towards the large harmonics when substruc- 
ture are accounted for. This effect is more pronounced when 
substructures have a cusp profile. 

Several improvements of the present investigation are 
envisioned since the ultimate goal of the method is to pro- 
vide a clear estimate of the amount of substructures in ob- 
servations. 

• Statistically, the properties of the cross-power spectrum 
of d/o(f?)/df? and fi(6) will be instructive. A clear charac- 
terisation of the covariance of these fields observed in Fig. 
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[3] along with the result of Alard (2008) upon which, as high 
multipole order n, about the same power is contained in 
each of the fields, would allow us to reduce the dimension- 
ality of the problem and perhaps only consider either fi(8) 
or df (6)/de . 

• In this vein, it remains to be shown to which extent the 
statistical properties of fi(6) and d/o(#)/d# can be ap- 
proximated as Gaussian random fields. If so, the realization 
of mock giant arcs would be greatly simplified. Concerning 
this point, the number of experiments of the present work 
must be increased to provide a clear diagnostic. 

• A natural extension of this work would be to consider 
more realistic lenses while taking into account the dynam- 
ics of substructures inside the host halo, its connection to 
cosmology via the expected statistical distribution of sub- 
structures (see e.g. Pichon & Aubert 2006), as well as star 
formation mechanisms. Indeed, stripping proccesses caused 
by tidal forces may lead to more complex structures. For in- 
stance, the study of some merger events in the phase-space 
(radial velocity versus radial distance) reveals the formation 
of structures quite similar to caustics generated in secondary 
infall models of halo formation (Peirani & de Freitas Pacheco 
2007). On the other hand, cooling processes and subsequent 
star formation may lead to steeper dark matter profile due 
to adiabatic contraction. Then, such processes should have 
an impact in the amplitude of high orders of our study. 

• Obviously a perturbation of the simple circularly sym- 
metrical case needs not lie in the same lens plane as the 
main deflecting halo. Uncorrelated halos superposed along 
the line of sight to the background source may well introduce 
perturbations on top of substructures belonging to the main 
halo. This will contribute to some additional shot noise back- 
ground in the power spectra of fi(0) and d/o(r?)/d# , which 
has to be quantified and subtracted off using ray-tracing 
through large simulated volumes. This work is beyond the 
scope of the present analysis. 

• On the path to a possible inversion yielding f\ (6) and 
dfo(6)/dO from observed arcs shapes and locations, several 
unknowns left on the rhs of Eq. (|33[l have to be controlled 
and will need to be fitted for in a non-linear way in order 
to attempt a reconstruction of fields fi(9) and d/o((9)/d(9 
. In addition, we only have considered a simple representa- 
tion of the background source. A lumpier background source 
will translate into a less regular arc with some small scale 
signature in the observable quantities such as dr±. How- 



ever, the replication of these internal fluctuations along the 
arcs and, possibly, in the counter image, as well as the in- 
formation contained in the various isophotes could allow us 
to reconstruct fi(0) and dfo(8)/d9 directly from the ob- 
servations. In this respect, Alard (2008) provides a general 
inversion method when two circular sources are considered 
for instance. 

To conclude, the upcoming generation of high spatial 
resolution instruments dedicated to cosmology (e.g. JWST, 
DUNE, SNAP, ALMA) will provide us with an unprece- 
dented number of giant arcs at all scales. The large samples 
expected will make standard lens modellings untractable 
and require the development of new methods able to capture 
the most relevant source of constraints for cosmology. In this 
respect, the perturbative method we present here may turn 
out to be a promising research line. 



6 ACKNOWLEDGEMENTS 

S. P. acknowledges the financial support through a 
ANR grant. This work made use of the resources avail- 
able within the framework of the horizon collaboration: 
http://www.projet-horizon.fr. It is a pleasure to thanks 
T. Sousbie, K. Benabed, S. Colombi, B. Fort and G. Lavaux 
for interesting conversations. We thank the referee for his 
useful comments that helped to improve the text of this 
paper. We would also like to thank D. Munro for freely 
distributing his Yorick programming language (available at 
http://yorick.sourceforge.net/) which was used during 
the course of this work. 



REFERENCES 

Amara A., Metcalf R. B., Cox T. J., Ostriker J. P., 2006, 

MNRAS, 367, 1367 

Alard, C. 2007, MNRAS, 382, L58 

Alard, C. 2008, MNRAS, 388, 375 

Aubert, D., Amara, A., & Metcalf, R. B. 2007, MNRAS, 
376, 113 

Aubert, D., Pichon, C, & Colombi, S. 2004, MNRAS, 352, 
376 

Bartelmann, M., Steinmetz, M., & Weiss, A. 1995, A&A, 
297, 1 



14 S. Peirani et al. 



Belokurov, V., et al. 2007, ApJ, 654, 897 
Bertschinger, E. 2001, ApJS, 137, 1 

Bode, P., Ostriker, J. P., & Turok, N. 2001, ApJ, 556, 93 
Bradac M., Schneider P., Steinmetz M., Lombardi M., King 
L. J., Porcas R., 2002, A&A, 388, 373 
Bradac M., Schneider P., Lombardi M., Steinmetz M., 
Koopmans L. V. E., Navarro J. F., 2004, A&A, 423, 797 
Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2001, 
ApJ, 548, 33 

Bullock, J. S., Kolatt, T. S., Sigad et al, 2001, MNRAS, 
321, 559 

Cimatti, A., Daddi, E., & Renzini, A. 2006, A&A, 453, L29 
Chen J., Rozo E., Dalai N., Taylor J. E., 2007, ApJ, 659, 
52 

Cole, S., & Lacey, C. 1996, MNRAS, 281, 716 

Cole, S., et al. 2005, MNRAS, 362, 505 

Coh'n, P., Klypin, A. A., Kravtsov, A. V., & Khokhlov, 

A. M. 1999, ApJ, 523, 32 

Coh'n, P., Avila-Reese, V., & Valenzuela, O. 2000, ApJ, 
542, 622 

Croft, R. A. C, Weinberg, D. H., Bolte, M., Buries, S., 
Hernquist, L., Katz, N., Kirkman, D., & Tytler, D. 2002, 
ApJ, 581, 20 

Dalai N., Kochanek C. S., 2002, ApJ, 572, 25 

de Blok W. J. C, McGaugh S. S., Bosma A. & Rubin V. 

C, 2001, ApJ 552, L23 

De Lucia, C, Kauffmann, C, Springel, V., White, S. D. M., 
Lanzoni, B., Stoehr, F., Tormen, C, & Yoshida, N. 2004, 
MNRAS, 348, 333 

Diego, J.M.,Protopapas, P., Sandvik, H.B., Tegmark, M., 
2005, MNRAS, 360, 477 

Dolag, K., Bartelmann, M., Perrotta et al., 2004, A&A, 
416, 853 

Fedeli, C, Meneghetti, M., Bartelmann, M., Dolag, K., & 
Moscardini, L. 2006, A&A, 447, 419 

Gao, L., De Lucia, C, White, S. D. M., & Jenkins, A. 2004, 
MNRAS, 352, LI 

Gavazzi, R., Treu, T., Rhodes, J. D., Koopmans, L. V. E., 
Bolton, A. S., Buries, S., Massey, R. J., & Moustakas, L. A. 
2007, ApJ, 667, 176 

Gentile, G., Salucci, P., Klein, U., Vergani, D., & Kalberla, 
P. 2004, MNRAS, 351, 903 

Ghigna, S., Moore, B., Governato, F., Lake, G., Quinn, T., 
& Stadel, J. 1998, MNRAS, 300, 146 

Ghigna, S., Moore, B., Governato, F., Lake, G., Quinn, T., 
& Stadel, J. 2000, ApJ, 544, 616 

Glazebrook K. et al, 2004, Nat 430, 181, |astro-ph/0401037| 
Hayashi, E., Navarro, J. F., Taylor, J. E., Stadel, J., & 
Quinn, T. 2003, Apj, 584, 541 

Hennawi, J. F., Dalai, N., Bode, P., & Ostriker, J. P. 2007, 
ApJ, 654, 714 

Hernquist, L. 1990, ApJ, 356, 359 

Horesh, A., Ofek, E. O., Maoz, D., Bartelmann, M., 
Meneghetti, M., & Rix, H.-W. 2005, ApJ, 633, 768 
Kasun, S. F., & Evrard, A. E. 2005, ApJ, 629, 781 
Kauffmann G., White S. D. M. & Guiderdoni B., 1993, 
MNRAS 264, 201 

Kazantzidis, S., Mayer, L., Mastropietro, C, Diemand, J., 

Stadel, J., & Moore, B. 2004, ApJ, 608, 663 

Kazantzidis, S., Kravtsov, A. V., Zentner, A. R., Allgood, 

B. , Nagai, D., & Moore, B. 2004, ApJL, 611, L73 
Keeton C. R., Gaudi B. S., Petters A. O., 2003, ApJ, 598, 



138 

Klypin A., Kravtsov A. V., Valenzuela O. & Prada F., 1999, 
ApJ 522, 82 

Knebe, A., Draganova, N., Power, C, Yepes, G., Hoff- 
man, Y., Gottloeber, S., & Gibson, B. K. 2008, astro- 
ph/0802.1917 

Kochanek, C. S., Keeton, C. R., & McLeod, B. A. 2001, 
ApJ, 547, 50 

Kochanek C. S., Dalai N., 2004, ApJ, 610, 69 
Koopmans L. V. E., 2005, MNRAS, 363, 1136 
Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., Klypin, 
A. A., Gottlober, S., Allgood, B., & Primack, J. R. 2004, 
ApJ, 609, 35 

Leauthaud, A., et al. 2007, ApJS, 172, 219 
Li, G.-L., Mao, S., Jing, Y. P., Bartelmann, M., Kang, X., 
& Meneghetti, M. 2005, ApJ, 635, 795 
Maccio A. V., Moore B., Stadel J., Diemand J., 2006, MN- 
RAS, 366, 1529 

Mao S., Jing Y., Ostriker J. P., Weller J., 2004, ApJ, 604, 
L5 

Meneghetti M., Argazzi R., Pace F., Moscardini L., Dolag 
K., Bartelmann M., Li G., Oguri M., 2007a, A&A, 461, 25 
Meneghetti M., Bartelmann M., Jenkins A., Frenk C, 
2007b, MNRAS, 381, 171 

Moore B., Ghigna S., Governato F., Lake G., Quinn T., 
Stadel J. & Tozzi P., 1999, ApJ 524, L19 
Nagai, D., & Kravtsov, A. V. 2005, ApJ, 618, 557 
Natarajan P., De Lucia G., Springel V., 2007, MNRAS, 
376, 180 

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 
462, 563 

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 
490, 493 

Oguri, M., Lee, J., & Suto, Y. 2003, ApJ, 599, 7 
Palunas P. & Williams T.B., 2000, AJ 120, 2884 
Peirani, S., & de Freitas Pacheco, J. A. 2007, 
|astro-ph/0701292| 

Pereira, M. J., Bryan, G. L., & Gill, S. P. D. 2008, ApJ, 
672, 825 

Pichon, C, & Aubert, D. 2006, 368, 1657 

Prunet, S., Pichon, C, Aubert, D., Pogosyan, D., Teyssier, 

R., & Gottloeber, S. 2008, astro-ph/0804.3536 

Salucci, P., & Burkert, A. 2000, ApJl, 537, L9 

Shaw, L. D., Weller, J., Ostriker, J. P., & Bode, P. 2006, 

ApJ, 646, 815 

Spergel, D. N., & Steinhardt, P. J. 2000, Physical Review 
Letters, 84, 3760 

Springel, V., White, S. D. M., Tormen, G., & Kauffmann, 
G. 2001, MNRAS, 328, 726 

Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 
361, 776 

Springel, V. 2005, MNRAS, 364, 1105 

Stoehr, F., White, S. D. M., Springel, V., Tormen, G., & 

Yoshida, N. 2003, MNRAS, 345, 1313 

Suyu S. H., Blandford R. D., 2006, MNRAS, 366, 39 

Vale A., Ostriker J. P., 2004, MNRAS, 353, 189 

van den Bosch, F. C, et al. 2007, MNRAS, 376, 841 

Viel M., Matarrese S., Theuns T., Munshi D. & Wang Y., 

2003, MNRAS 340, L47 

Viel M., Haehnelt M.G. & Springel V., 2004, MNRAS 354, 
684 

Warren, M. S., Quinn, P. J., Salmon, J. K., & Zurek, W. H. 



Numerical investigation of lenses with substructures using the perturbative method. 



15 



1992, ApJ, 399, 405 

Yoo, J., Kochanek, C. S., Falco, E. E., & McLeod, B. A. 

2005, ApJ, 626, 51 

Yoo, J., Kochanek, C. S., Falco, E. E., & McLeod, B. A. 

2006, ApJ, 642, 22 



